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ABSTRACT 


The  research  described  in  this  report  is  aimed  at 
characterizing  the  non-equilibrium  flow  and  coupled 
chemistry  associated  with  the  interaction  of  rocket  exhaust 
effluents  in  the  upper  atmosphere.  Several  aspects  of  the 
current  study  should  be  noted.  Although  continuum  fluid 
descriptions  based  on  the  Navier-Stokes  equations  are 
valid  at  low  altitudes,  at  high  altitudes  rarified  effects 
must  be  considered.  To  investigate  this  transition  and 
the  ultimate  analytical  modeling  of  rarified  flows,  a 
Monte  Carlo  technique  applied  to  a  Boltzmann  equation 
formulation  is  being  studied.  The  technique  is  described 
in  this  report,  together  with  the  progress  to  date. 

The  major  observable  from  rocket  exhaust-atmospheric 
interactions  is  radiation  in  various  wavelength  regions. 

To  understand  the  origin  of  this  emission  we  have  focused 
attention  on  the  emissions  as  generated  from  the  exothermicity 
of  chemical  reactions  involving  unburned  fuel  species  and 
0  and  OH  radicals.  These  experimental  studies  have  shown 
significant  emissions,  particularly  in  the  infrared,  as  a 
result  of  direct  product  vibrational  excitation  of  such 
species  as  OH,  HjO  and  CO2.  The  specific  measurements  and 
results  are  described  in  Section  II  of  this  report. 

Lastly,  in  order  to  characterize  the  non-equilibrium 
chemistry  in  expanding  gas  flow  fields,  a  detailed  knowledge 
of  relaxation  processes  must  be  obtained.  Together  with 
available  information,  both  theoretical  calculations  and 
laboratory  measurements  are  herein  described  to  characterize 
the  relaxation  processes  pertinent  to  expanding  combustion 
products  in  a  rarified  atmosphere. 


re  a 


PREFACE 


The  research  described  in  this  report  falls  into 
four  main  categories/  identified  in  the  contractual 
work  statement  as  follows: 

a)  Calculate  the  continuum  flow  structure  and  rate  chemistry 
of  high  altitude  plumes. 

b)  Calculate  the  vibrational  excitation  cross-sections  for 
plume  relevant  species  using  available  interaction  potential 

data. 

c)  Measure  the  quantum  yields  in  the  infrared  region  for  the 
chemi-excitation  involving  0  or  OH  and  plume  species. 

d)  Measure  the  inter-molecular  potential  between  atmic  oxygen 
and  various  plume  species  as  an  input  to  the  calculation 
of  vibrational  excitation  cross  sections. 

Items  (a)  and  (b)  are  discussed  in  Section  I,  item  (c)  in 
Section  II,  and  item  (d)  in  Section  III. 
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SECTION  I 


Fluid  dynamics 

In  this  portion  of  the  research  we  are  concerned  with 

developing  and  testing  procedures  for  computing  the  flow 

fields  of  high-altitude  rocket  exhaust  plumes  and  of  atmospheric 

releases  which  may  be  of  interest  in  connection  with  simulations 

of  plume  effects.  The  major  portion  of  an  exhaust  plume  from 

5 

a  typical  large  rocket  (thrust  about  10  lb)  is  in  a  continuum 

flow  regime  at  altitudes  below  about  160  km.  At  altitudes 

below  about  110-120  km. ,  the  mixing  of  the  exhaust  gases  with 

air  is  turbulent.  At  higher  altitudes  a  reverse  transition 

occurs,  and  the  plume  becomes  essentially  laminar,  although 

developing  instabilities  may  persist  for  a  few  seconds.  At 

these  altitudes  the  Navier-Stokes  equations  are  presumed  to 

be  sufficient  to  describe  local  conditions  in  the  exhaust 

plume  flow  field.  Procedures  for  solving  these  equations  for 

thrusting  plumes  with  zero  angle  of  attack  are  available  as 

1  2 

the  result  of  previous  investigations  '  and  exemplificative 
calculations  are  available  .  Under  the  current  contract,  work 
has  been  directed  at  examining  the  deviations  from  this 
description  at  altitudes  whore  rarefaction  effects  become 

important,  and  also  at  developing  some  comparisons  with 

4  .  . 

recent  experimental  data  .  The  procedures  for  examining 

rarefaction  effects  have  been  formulated  and  coded  for  computer 

solution,  but  have  as  yet  produced  no  significant  results. 

The  comparison  with  experiment  has  just  been  initiated. 

To  describe  rarefaction  effects,  we  are  making  use  of  an 
ingenious  and  physically  appealing  approach  to  numerical 

5 

solution  of  the  Boltzmann  equation  developed  by  Bird  .  The 
approach  used  is  reminiscent  of  the  molecular  dynamics  studies 
of  Alder  and  Wainwright  ,  in  that  the  solution  is  developed  as 
the  solution  to  a  Liouville  equation  for  a  small  sample  of 
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test  particles  which  represent  the  actual  particles  in  the 
real  situation.  Alder  and  Wainwright  evaluated  the  macroscopic 
properties  of  their  systems  by  computing  in  detail  the 
microscopic  motions  of  an  assembly  of  molecules  with  prescribed 
intermolecular  forces.  Macroscopic  properties  were  evaluated 
by  averaging  over  a  number  of  realizations  of  the  microscopic 
system.  The  power  of  available  computers  (this  work  was 
performed  in  the  1950's)  and  the  necessity  for  following  the 
motion  in  detail,  even  during  collisions  where  rapid  changes 
in  velocity  occur,  limited  the  scope  of  Alder  and  Wainwright' s 
work  primarily  to  consideration  of  phase  transitions  in 
systems  of  two  space  dimensions,  and  to  a  small  sample  size. 

Bird's  approach  is  similar  to  that  of  Alder  and  Wainwright, 
in  that  the  motion  of  a  set  of  test  particles  is  followed  in 
detail  between  collisions,  but  differs  in  that  collisions  are 
not  evaluated  in  detail,  instead.  Bird  assumes  that  collisions 
are  very  short-lived  on  the  time  scale  of  interest  for  the 
overall  solution,  and  evaluates  their  effects  by  a  random 
sampling  procedure  or  Monte  Carlo  technique.  (He  calls  his 
procedure  the  "Monte  Carlo  direct  simulation  method.")  The 
procedure  is  best  illustrated  by  an  example. 

Suppose  that  we  wish  to  calculate  conditions  in  a  rarefied 
gas  in  a  region  lying  within  a  bounding  surface  at  which  the 
state  of  the  gas  (i.e.,  the  local  density  and  velocity  distribution 
function)  is  prescribed.  The  region  within  the  bounding  surface 
is  divided  into  a  number  of  cells  whose  dimensions  are  of  the 
order  of  a  mean  free  path.  Initial  conditions  are  prescribed 
in  each  cell,  and  each  is  assigned  a  set  of  particles  repre¬ 
senting  the  gas  molecules.  The  number  of  these  particles  in 
each  cell  depends  upon  the  assigned  density.  Their  positions 
within  the  cell  are  assigned  at  random,  and  their  velocities 
are  chosen  by  a  random  sampling  procedure  consistent  with  the 
initially  assigned  velocity  distribution.  The  motion  of  these 
particles  is  then  followed  over  a  succession  of  time  steps. 
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At  the  beginning  of  each  time  step,  changes  in  the 
velocities  of  the  individual  particles  due  to  collisions  are 
evaluated.  A  particle  is  selected  at  random  from  those  in  a 
given  cell.  A  second  particle  is  chosen  at  random  from  the  same 
cell  as  the  first.  This  pair  of  particles  is  either  assumed 
to  collide  or  assumed  not  to  collide;  the  probability  of 
collision  is  proportional  to  the  relative  velocity  of  the 
collision  partners.  If  the  pair  are  to  collide,  a  set  of 
collision  parameters  is  chosen  at  random  and  the  resulting 
new  velocities  after  this  collision  are  evaluated.  The 
colliding  partners  are  assigned  these  new  velocities.  A  new 
set  of  collision  partners  is  then  chosen  in  that  cell,  and  the 
calculation  of  collisions  is  carried  on  in  this  fashion  until 
a  number  appropriate  to  the  local  collision  frequency  and  the 
time  step  have  occurred.  One  then  proceeds  to  the  next  cell 
and  repeats  this  process.  When  an  appropriate  number  of 
collisions  in  each  cell  have  been  evaluated,  all  the  particles 
are  moved  through  a  distance  given  by  the  product  of  their 
velocities  and  the  time  step.  If  a  particle  crosses  a  cell 
boundary,  it  is  removed  from  the  calculation.  New  particles 
enter  the  calculation  at  random  from  points  on  the  boundary  at 
an  average  rate  and  with  properties  determined  from  the  local 
boundary  conditions.  The  calculation  of  collisions,  new 
positions,  particle  losses,  and  particle  gains  is  repeated  at 
each  time  step. 

Periodically,  the  calculation  is  sampled.  The  density  in 
a  cell  is  proportional  to  the  number  of  test  particles  it 
contains  and  inversely  proportional  to  its  volume.  The  velocity 
distribution  function  is  developed  as  a  histogram  of  the 
velocities  of  the  test  particles.  No  information  is  requested 
on  a  scale  smaller  than  the  cell  size. 
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Two  kinds  of  problems  can  be  treated  with  this  method. 

Steady  rarefied  flows  can  be  computed  by  imposing  constant 
boundary  conditions  and  allowing  the  steady- state  solution 
to  evolve  as  the  long-time  limit  of  a  time-dependent  calculation. 
In  any  particular  sample,  the  deviation  of  the  density  and 
velocity  distribution  from  their  average  values  may  be  large; 
more  correct  values  are  obtained  by  time-averaging  these 
properties  in  each  cell  after  a  (statistically)  steady  state 
has  been  attained.  Unsteady  flows  are  treated  as  time-dependent 
calculations,  and  the  scatter  of  the  results  is  reduced  by 
averaging  over  a  number  of  successive  calculations  for  the 
same  problem. 

This  procedure  was  originally  applied  by  Bird  to  calculating 

flow  conditions  in  shock  waves  and  to  calculations  of  drag 

coefficients  in  the  transition  flow  regime.  It  has  been  used 

7  8 

by  Vogenit  z,  et  al.  '  for  calculating  drag  coefficients  on 
high-altitude  bodies  and  to  investigate  some  aspects  of  the 
transition  from  continuum  to  rarefied  flow  on  simple  shapes. 

The  procedure  can  be  regarded  as  an  algorithm  for  solving  the 
Boltzmann  equation.  The  arguments  for  regarding  it  as  such 
have  been  put  forward  by  Bird  in  a  recent  paper.  Basically, 
the  procedure  gives  one  all  the  information  which  one  asks 
from  the  Boltzmann  equation  to  an  accuracy  which  depends  upon 
the  sample  size,  including  the  number  of  test  particles  and  the 
length  of  the  interval  over  which  a  time  average  is  taken.  The 
random  sampling  of  collisions  is  consistent  with  the  level  of 
information  contained  in  solutions  to  the  Boltzmann  equation  so 
long  as  this  sampling  is  done  in  a  way  consistent  with  collision 
probabilities  in  the  real  flow  being  simulated. 

We  are  in  the  process  of  adapting  Bird's  procedure  to 
calculating  flows  in  high-altitude  atmospheric  releases  and 
(in  the  small-disturbances  approximation)  to  rocket  exhaust 
plumes.  The  procedure  will  ultimately  be  used  as  a  time-dependent 
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calculation  in  one  space  and  two  velocity  dimensions.  (At 
the  present  time,  the  applications  have  been  limited  to 
steady-state  flows.)  A  computer  code  has  been  prepared  to 
evaluate  the  Monte  Carlo  solution  of  rarefied  flow  problems. 

The  code  in  its  current  form  is  described  below.  Following 
this  description  we  discuss  planned  improvements  which  will 
allow  us  to  investigate  upper-atmosphere  releases  and  chemical 
reactions. 

We  initiate  the  calculation  by  specifying  cell  boundary 
positions  along  the  x  axis.  The  cell  dimensions  are  currently 
selected  to  make  them  of  the  order  of  the  expected  mean  free 
path  in  that  location  after  the  calculation  has  achieved  the 
statistically  steady  state.  Each  cell  is  assigned  an  initial 
temperature,  number  density,  and  mean  velocities  parallel  and 
transverse  to  the  x  axis.  The  calculation  assigns  to  each 
cell  a  number  of  particles  proportional  to  the  assumed  number 
density.  The  particles  are  located  at  random  between  the 
boundaries  of  the  cell.  The  code  assigns  these  positions  by 
calling  a  random  number  generator  subroutine  (RANDU  on  the  IBM  360) 
which  produces  a  random  number  between  0  and  1;  the  position  is 
thus 

xp  =  *i  +  R(x2  -  xx) 

where  x2  and  x^  are  the  upper  and  lower  boundaries  of  the  cell 

and  R  is  the  random  number.  The  initial  velocities  assigned 

are  selected  at  random  in  a  manner  consistent  with  a  Maxwellian 

distribution  at  the  local  temperature  and  with  the  mean 

velocities.  Thus  if  the  local  temperature  is  T  and  the  mean 

velocity  is  u  =  iux  +  ^uy  +  £uz,  the  velocities  of  the  particles 

may  be  assigned  by  using  the  IBM  subroutine  GAUSS,  which 

produces  random  numbers  with  a  Gaussian  distribution  about  a 

specified  mean  value  with  a  specified  standard  deviation  -  in 
2 

this  case,  a  =  kT/m. 
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The  temperatures  and  velocities  outside  of  the  first  and 
last  cells  in  the  calculation  are  assumed  constant.  These 
boundary  conditions  complete  the  specification  of  the  problem. 
The  location  of  the  boundaries  is  chosen  far  enough  from  th  2 
region  of  interest  so  that  their  effects  are  minimal.  (For 
example,  consider  a  normal  shock  wave.  Here  we  may  assign 
at  the  upstream  boundary  T^  o^,  and  u^  and  at  the  downstream 
boundary  T2,  p2,  and  u2  according  to  the  Rankine-Hugoniot 
relations.  If  the  boundaries  are  chosen  too  close  together, 
the  shock  wave  will  be  distorted.  As  the  boundaries  are  moved 
further  apart,  the  solution  approaches  more  closely  the 
solution  obtained  when  the  boundary  conditions  are  imposed 
only  at  great  distances  from  the  shock.  A  reasonable  compromise 
between  computation  time  and  accuracy  is  usually  necessary.) 

The  calculation  is  initiated  by  calculating  collisions  in 
the  first  cell.  A  "molecule"  is  chosen  at  random  from  those 
in  the  cell,  and  a  second  molecule  is  also  selected  at  random. 

The  probability  that  these  two  molecules  will  collide  is 
proportional  to  their  relative  velocity,  which  can  be  determined 
from  the  assigned  velocity  components.  Thus  if 

!vrl  >  v0R, 

where  v0  is  a  large  velocity  and  R  a  random  number,  the  pair 
is  allowed  to  collide.  If  not,  a  second  pair  is  chosen  and  the 
process  is  repeated.  When  a  pair  has  been  selected,  the 
collision  dynamics  are  calculated  for  a  randomly  selected 
impact  parameter. 

In  the  current  version  of  the  code,  a  hard-sphere  potential 
has  been  chosen  for  simplicity.  Since  the  scattering  of  hard 
spheres  is  isotropic,  the  orientation  of  the  new  relative 
velocity  is  a  random  function  of  6  and  cp  on  the  unit  sphere. 

New  velocity  vectors  relative  to  the  center  of  mass  for  each 
particle  can  be  calculated  in  terms  of  two  random  numbers,  R^  and 
R2,by  setting 
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and 


cp  =  2rr 


sin  0  =  4Rj  (1  - 

and  calculating  new  velocity  components 


v  =  v  ±  kv  sin 0  cos  cp 
x  cm,x  ^  r 


v  =  v  ±  kv  sin0  cos  cp 
y  cm,y  r 


=  v  „  ±  cos  0 
cm,  z  r 


where  the  +  sign  is  chosen  for  one  particle  and  the  -  sign  for 
the  other. 

At  each  collision  in  a  particular  cell  a  time  counter  tc 
is  advanced  by  the  amount 

6tc  =  2/dLTT  a*  n.vr) 

where  n.  is  the  number  density  in  the  ith  cell,  N.  is  the 

actual  number  of  particles  in  the  simulation  calcuation,  and 

nag  is  the  collision  cross-section.  The  calculation  of 

collisions  in  the  cell  is  continued  until  the  time  counter  tc 

exceeds  t  ,  a  pre-selected  value  which  is  small  compared  to  the 

mean  collision  time.  (In  practice,  we  examine  the  chances  of 

t  exceeding  t  before  each  collision.  If  t  +  5t  >  t  ,  we 
c  n  c  c  11 

allow  the  final  collision  if 


c 


and  disallow  it  otherwise.)  When  t  exceeds  tn,  the  calculation 
proceeds  to  the  next  cell. 
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After  all  collisions  have  been  calculated,  in  all  cells, 
each  particle  is  moved  through  a  distance 


If  the  new  value  of  xi  lies  within  another  cell,  the  particle 
is  removed  from  the  list  of  particles  contained  in  the  cell 
it  left  and  added  to  the  list  of  particles  contained  in  the 
cell  it  entered,  if  the  particle  has  crossed  one  of  the 
boundaries,  it  is  removed  from  the  calculation. 

New  particles  are  allowed  to  enter  the  calculation  from 

upstream  and  downstream,  using  a  procedure  developed  by 
9 

Haviland  .  The  number  flow  of  particles  into  the  region  of 
calculation  across  the  upstream  boundary  is 


r  r 

J0  dv  v  J_« 

°  XX 


dv 


'  i-»  dv  f.  (v  , v  , v  ) 

y  z  1  x  y  z 


where  fj^  is  the  velocity  distribution  function  imposed  at  the 
upstream  boundary,  in  our  calculation  the  velocity  distribution 
is  currently  a  Maxwellian  about  the  mean  velocity  u^,  so  that 


fdvv  fdv  r'av 

1  *■'  o  xxJ-cb  yj-ae 


n. 


dv 
os  y 


z  172  3" 

"  C1 


2  2U 
+v  +v 

_ Z _ 


2 

z 


where  c^  =  ^kT^/m)^,  and  n^  is  the  number  density.  The 
integrals  are  easily  evaluated  to  give 
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T,2/  2 


r  ,  u,  -i 

i=  nici|>*'  6  +  2^x{1  +  erf(^> 


Likewise,  the  number  flow  into  the  region  of  calculation 
across  the  downstream  boundary  is 


o  r  x 

n2  =  n2c2[^T  6 


„2/  2 

-U2/c2 


u2  U 

4"  U  +  (4}  }J 


Since  particles  enter  at  random  across  the  upstream  boundary, 

the  probability  that  N  particles  enter  during  a  time  t 

th  ^ 

across  the  j  boundary  is 

N  0  i  0 

£  (t  n  ) 1  -t  n 

P  .  =  L  ----3  -  e  n  3 

n » i  i_Q  j 


so  that  the  number  of  new  particles  to  be  added  can  be  calculated 

by  setting  B  .  equal  to  a  random  number  R. 
n»  j 

Velocities  are  assigned  to  the  new  particles  according  to 
the  distribution  function  at  the  boundary  at  which  they 
entered.  For  a  Maxwellian  distribution  about  a  mean  velocity 
at  the  upstream  boundary,  the  components  in  the  y  and  z 
directions  are  assigned  by  using  GAUSS.  The  value  of  ux  is 
determined  by  solving  the  equation 
r®  /v-U,  \  2 


r*  f-“i\ 

.u  ve-l -5-7 

f^W2 


,U,-u  \2  n/tT  U,  i  /U.-U  v  v 

,-(-y)  4^)} 

,'(U1/°1)  +V-^i  {  1  +  erf  (c^)  } 


fU,  -u 
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where  R  is  a  random  number  and  the  solution  can  be 
a  Newton-Raphson  iterative  procedure  starting  with 
guess 


o 


obtained  by 
the  initial 


At  the  downstream  boundary,  a  like  procedure  is  followed  with 
replaced  by  -l^  and  c^  by  c The  solution  is  then 
multiplied  by  -1  to  give  the  correct  direction.  The  position 
of  each  newly  added  particle  is  determined  by  letting 


x.  =  x .  +  R  U  .  t  . 
i  j  x, i  m 


The  above  procedure  has  been  coded  for  the  IBM  360  and 
is  currently  being  checked  out.  At  the  current  time,  we  are 
applying  the  procedure  to  the  calculation  of  a  steady,  one- 
dimensional  plane  shock  wave.  We  have  as  yet  been  unable  to 
achieve  convergence  to  the  final  result  (attainment  of  a 
statistically  steady  state) .  The  running  time  on  our  computer 
(IBM  360-75)  is  rather  long.  The  shock  calculation  comprises 
20  cells  and  about  500  particles,  and  at  one  hours'  running 
time  had  not  yet  relaxed  to  statistically  steady  densities, 
mean  velocities,  and  temperatures.  No  extensive  attempts  to 
optimize  the  code  have  yet  been  made,  although  it  is  clear 
that  to  obtain  a  calculation  which  results  in  convergence  of 
high  moments  of  the  distribution  function  (a  necessary  condition 
for  investigation  of  effects  on  chemical  reactions  of  a  non¬ 
continuum  environment)  such  optimization  will  be  needed.  We 
are  now  in  the  process  of  making  it  possible  to  write  out 
intermediate  output  on  tape,  so  that  a  calculation  may  be  run 
in  several  steps  of  moderate  individual  length  rather  than  in 
one  very  long  step. 


We  plan  to  employ  the  foregoing  procedure  to  assess  the 
effects  of  rarefaction  upon  reaction  and  excitation  processes 
in  hiqh-altitude  plumes.  The  calculation  (in  cylindrical 
geometry)  is  actually  a  small-disturbance  approximation  to  the 
plume  flow  field;  in  comparing  results  with  the  plume  calculation, 
this  must  be  taken  into  account.  In  order  to  be  able  to 
assess  the  effects  on  reactions  independently  of  the  effects  of 
the  small-disturbance  calculation,  we  felt  it  desireable  to 
develop  a  small-disturbance  code  for  the  continuum  flow.  This 
code  should  also  be  useful  for  calculating  conditions  in  high- 
altitude  releases.  The  code  is  based  upon  a  Lagrangian  finite- 
difference  procedure  very  similar  to  that  employed  in  MULTITUBE, 
and  essentially  is  an  adaptation  of  MULTITUBE  to  a  one¬ 
dimensional  unsteady  flow. 
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The  equations  of  overall  continuity,  momentum ,  energy 
and  species  conservation  are 


+  P  (V  •  v)  =  0 


_  uv 

p  at 


=  -  VP  -  (V-t) 


=  -  V  •  q  -  V 


Z  hi  3i  - 


P  V -V  -  t  :  v 


p  5t~ 


*  -V.J.  +  w. 


J 

Here  the  time  derivative  ^  is  taken  following  a  fluid  element. 
The  heat  flux  q  is  -kVT,  and  other  quantities  are  defined  in 
the  table  of  nomenclature. 

We  now  specialize  to  a  system  in  one  space  dimension  r , 
and  introduce  a  metric  coefficient  j ,  which  takes  on  the  values 
0,1,  and  2  for  planar,  cylindrical,  or  spherical  symmetry,  re¬ 
spectively.  All  fluxes  and  velocities  are  zero  except  in  the 
radial  direction  and  we  let  v  =  vr,  q  =  qr,  and  The 

non-vanishing  parts  of  the  stress  tensor  are  the  diagonal  com¬ 


ponents 


“  -p[2  i  - 17'5] 
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Then 


and 


T  :7v 


3v  . 
Trr  37  +  3 


v 

r  Tee  * 


We  can  then  write  the  governing  equations  in  terms  of  a  "radial" 
normal  stress  Trr  and  an  "angular"  normal  stress  tQq  as 


dp 

3t 


^  h  (rjv) 


dv 

p  at 


dE 

p  at 


_  _  3P  _  1  3  (rjT  \  _  i  T 
"  37  IT  3r  l  n)  r  T 


00 


=  ‘?Mr3(q+£HiJi)| 


(rJv) 


9v  _  jv  T 
Trr  3r  “  r  00 


dc . 
p  3F" 


r  3r  l  iJ  i 


We  will  develop  the  Lagrangian  finite-difference  approximation 
to  these  equations  in  tarms  of  a  model  consisting  of  an  array 
of  expandable  cells,  with  movable  boundaries  separating  them. 

It  is  natural  to  define  some  quantities  as  cell  quantities  and 
others  as  boundary  quantities.  The  Lagrangian  velocity  v  is 
naturally  associated  with  a  moving  boundary,  for  example,  while 
the  pressure  and  the  internal  energy  are  naturally  associated 
with  a  cell.  We  can  then  use  the  momentum  equation  to  describe 
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the  motion  of  the  cell  boundaries,  and  the  remaining  equations 
to  describe  the  properties  of  the  fluid  in  the  cells.  Formally, 
we  proceed  by  introducing  a  coordinate  transformation 

dy  =  c jr^  dr  , 


where  y  has  the  dimensions  of  a  (j+1) -dimensional  volume  and 
Cj  is  1,  2tt,  or  4tt  as  j  is  0,  1,  or  2.  Then 


3_ 

3r 


V 


3_ 

3y 


and  the  equations  become 
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We  now  pass  from  this  transformed  set  of  equations  to  the 
finite-difference  equations.  Let  the  flow  be  divided  into  a  num¬ 
ber  of  cells,  and  define  the  boundary  of  each  cell  in  such  a  way 
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that  the  total  mass  within  each  cell  does  not  change  with  time. 

The  cell  boundaries  are  then  identifiable  as  mean  particle  tra¬ 
jectories  in  the  moving  fluid.  While  individual  species  may 
pass  across  the  mean  trajectory  by  diffusion,  no  net  transfer  of 
mass  occurs.  The  mean  trajectories  are  analogous  to  streamlines 
in  flows  with  two  space  dimensions.  It  is  natural  to  define  the 
fluid  velocity  at  the  cell  boundaries,  and  thermodynamic  state 
properties  (temperature,  pressure,  composition,  internal  energy, 
density)  within  the  cells.  The  heat  and  mass  fluxes  are  also 
conveniently  defined  at  the  boundaries.  The  stresses,  which  are 
in  effect  additional  pressure  terms,  are  more  conveniently  de¬ 
fined  in  the  cells.  We  introduces  two  difference  operators,  Ak, 
which  takes  the  difference  in  boundary  properties  across  the  k 
cell,  and  $k,  which  takes  the  difference  in  cell  quantities  across 
the  kth  boundary.  We  then  adopt  the  convention  that  the  k  boun¬ 
dary  bounds  cell  k  on  the  "outside,"  i.e.,  the  direction  of  in¬ 
creasing  r.  Then  if  we  set  the  mass  of  a  cell  =  Pk^k 
the  equations  for  the  rates  of  change  of  the  cell  quantities 
become 


^  = 


constant 
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The 


m  dCi'k 
mk  cTt 


=  -e.A.(r^J.)  +  w.  ~ 

3  k\  1/  i  pk 


where  ^  and  ( ^  )• 

momentum  equation  and  the  definition  of  the  velocity  describe  the 
motion  of  the  cell  boundaries : 


= -dvk 
%  dt 


=  -c . 
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r,.3«,.P  +  6k(x\jJ 


‘k  ukJ 


+  2-7: 


rk  ee,k 


where  m^.  2  (™k  +  ^Yn)  and  T  qq  "  2(T60,k  +  T00,k+l)  * 

The  above  equations  provide  the  changes  in  cell  boundary  position 
and  velocity,  and  the  changes  in  composition  and  internal  energy 
within  a  cell,  over  a  small  time  increment.  The  thermal  and 
caloric  equations  of  state 


and 


E  =  E(T,p) 
P  =  P(T,p) 


provide  the  remaining  information  required  for  a  complete  speci¬ 
fication  of  state  properties  in  the  expanding  flow.  In  general 
we  will  consider  ideal  mixtures  of  perfect  gases,  so  that  these 
equations  become 
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and 

pRT  = 

A  solution  to  this  set  of  equations  consists  of  the  positions 

and  velocities  of  the  cell  boundaries,  together  with  the  state 
properties  in  the  cells,  as  functions  of  time  following  some 
prescribed  initial  condition.  The  procedure  for  evaluating  the 
solution  is  the  following:  Using  the  initial  velocities  and 
state  properties,  the  stresses,  the  heat  and  mass  fluxes,  and 
the  reaction  rates  within  the  cells  are  evaluated.  From  the 
calculated  stresses  and  the  known  pressures,  the  accelerations 
of  the  cell  boundaries  are  calculated,  using  the  equations  of 
motion.  After  a  small  time  step  6t,  new  positions  of  the  cell 
boundaries  are  computed  and,  from  the  new  cell  volumes  and  the 
constant  cell  masses,  new  cell  densities  are  evaluated.  The 
changes  in  composition  within  the  cells  is  then  evaluated  from 
the  previously  calculated  mass  fluxes  and  reaction  rates ,  and 
the  internal  energies  within  the  cells  are  computed  using  the 
energy  equation.  Finally,  from  the  now  known  density,  internal 
energy,  and  composition,  the  state  properties  are  determined 
using  the  equations  of  state.  The  flow  field  at  the  end  of  a  time 
step  <St  is  now  completely  specified,  and  the  calculation  is  re¬ 
peated  at  the  next  time  step.  The  calculation  marches  in  time 
in  this  fashion  until  the  entire  flow  has  been  evaluated. 
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One  can  show  by  summation  over  all  the  cells  that  the  difference 
scheme  identically  satisfies  the  integral  conservation  equations 
of  mass  and  energy. 

The  foregoing  procedure  has  been  coded  and  initial 
computer  calculations  are  now  being  performed  for  simple  free 
expansions  with  cylindrical  and  spherical  symmetry.  In 
preparing  the  code,  we  have  made  considerable  use  of  routines 
from  MULTITUBE.  Because  of  the  considerable  similarity 
between  these  problems,  most  of  the  FORTRAN  statements  from 
MULTI TUBE  could  be  used  directly  in  the  new  1-D  unsteady  code 
with  some  redefinition  of  the  variables.  Some  reprogramming 
was  necessary  in  order  to  account  for  differences  in  handling 
of  shocks  and  stresses.  Since  the  new  code  is  still  being 
checked  out,  we  shall  not  present  a  complete  description  here; 
we  expect  to  document  the  code  at  a  future  date. 

The  third  effort  under  the  fluid  mechanics  effort  is  the 
comparison  of  exhaust  plume  calculations  made  with  MULTITUBE 
to  experiments4  recently  performed  at  AEDC.  Specifically,  a 
calculation  is:  currently  under  way  for  the  test  conditions  IC5 
of  Reference  4,  which  consisted  of  a  helium  plume  for  a  1.61:1 
area  ratio  nozzle  injected  into  an  N2  free  station  at  Mach  18.15. 
These  calculations  and  the  comparison  with  experiment  (impact 
pressure  and  species  density)  will  be  the  subject  of  a  separate 
report. 
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SECTION  II 


I  ntroduction/Summary 

Prior  to  quantitative  model  predictions  on  the  radiative 
properties  of  exhaust  products  in  the  upper  atmosphere, 
several  important  chemical  reaction  rate  coefficients  must 
be  determined.  In  particular,  our  studies  have  focussed  on 
simple  exothermic  oxidation  reactions  involving  0  and  OH 
attack  on  various  hydrocarbon  fuel  species  and  accompanying 
vibrational  relaxation  processes. 

The  interest  in  simple  exothermic  reactions  is  directed 

toward  the  release  of  the  available  energy  of  the  reaction 

into  radiative  modes  of  the  product  molecules.  Measurements 

o 

in  our  laboratory  have  been  made  from  2000  A  to  25p,  in  an 
attempt  to  determine  the  reaction  mechanism  and  available 
radiative  modes. 

OH  is  a  strong  emitter  at  3064  A  (OH  A2  -*  X2tt)  as 
well  as  in  the  Meinel  Bands  from  7000  to  9000  A  and  in  the 
near  IR  at  2.7p,.  OH  emission  has  been  observed  as  a  result 
of  0  atom  attack  on  virtually  all  hydrocarbons  tested  to  date. 
In  addition,  we  find  that  the  subsequent  OH  attack  on  CH20 
results  in  strong  H20  emission  at  2.7y  as  well  as  some  6.3y 
radiation. 

Our  early  results  have  focused  on  the  identification 
and  determination  of  concentration  values  through  the  visible 
chemiluminescense  associated  with  the  attack  of  0  atoms  and 
OH  on  reactive  hydrocarbons.  Preliminary  survey  measurements 
have,  however,  been  made  out  to  25y  and  the  results  from  these 
measurements  are  discussed  herein.  As  the  mass  spectrometer 
system  is  completed  and  fully  incorporated  into  the  flow 
system,  the  current  stress  on  visible  emission  will  no  longer 
be  necessary. 
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Experimental 


All  studies  were  performed  in  the  discharge  flow  system 
shown  schematically  in  Figure  1.  Molecular  oxygen  (American 
Cryogenics,  purity  99.6%)  was  discharged  in  an  Evenson  or 
Broida  Cavity  at  2450  MHZ  using  a  Raytheon  PGM  10x2  diathermy 
system.  The  discharge  was  pulsed  at  238  Hz  for  simultaneous 
infrared,  visible  and  ultraviolet  measurements.  The  discharged 
oxygen  was  passed  into  the  flow  tube  reactor  through  a 
teflon  tube  with  a  radially-multiperforated  exit  stream  which 
permits  excellent  mixing  with  a  second  gas  and  prevents 
discharge  region  radiation  from  reaching  the  observation 
ports.  The  stainless  steel  reactor  is  lined  with  a  teflon 
sleeve  which  largely  prevents  wall  recombination  of  the  atomic 
species.  The  inside  diameter  of  the  reactor  is  one  inch  and 
typical  gas  velocities  are  30  to  40  m/sec  through  a  Heraeus 
Engelhard  E  -  225  mechanical  pump.  Six  ports  permit  radial 
viewing  through  KBr  windows  using  an  RCA  C31025C  uncooled 
photomultiplier  on  a  Jarrell  Ash  scanning  quarter  meter 
monochrometer  with  a  Keithley  602  Electrometer  recorded  on 
a  Honeywell  Electronik  194  recorder.  Signal  to  noise  ratios 
were  such  that  no  damping  or  time  integration  of  signals 
was  necessary  in  the  u.v.  and  visible.  Measurements  in  the 
infrared  were  made  using  a  Santa  Barbara  uncooled  PbS 
detector  (to  3 y,)  ,  Ge: Au  detector  (to  8y.)  and  Ge:Cu  detector 
(to  25y.) .  All  infrared  measurements  were  made  with  a  PAR  HR8 
phase  Sensitive  Lock  In  Amplifier1 . 

The  hydrocarbons  employed  as  reactant  gases  were  methane, 
ethane,  propane,  ethylene,  propylene,  trans-butene,  acetylene 
and  formaldehyde.  Linde  research  grade  methane  and  ethane 
were  employed.  Union  Carbide  chemically  pure  propane  was 
used.  The  propylene  and  trans-butene  from  Cryogenic  Sales 
had  purities  of  95%  and  99%  respectively.  Several  runs  were 
made  with  "house"  methane  with  no  perceptibly  different  results 
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Figure  1.  Schematic  of  discharge  flow  tube  reactor  and  detection 
system.  Legend: 

M  *  Monochromator  -  Jarrell  Ash,  quarter  meter, 
scanning. 

P  =.  Photomultiplier  -  RCA  C31025C  (uncooled) 

E  *  Electrometer,  D.C.  Amplifier  -  Keithley  602 

Pi  *  Pi  coampere  source  -  Keithley  261  , 

R  *  Recorder  -  Honeywell  Electronik  194 

Hi  =  High  voltage  D.C.  power  supply  -  Power  Designs 
Pacific  model  2K-10 
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Fig.  1-Schematic  of  Discharge  Flow  Tube  Reactor 

and  Detection  System 


Oxygen  atom  concentrations  were  determined  through  gas 
phase  NO2  titration;  typically  1%  atoms  were  present.  The 
spectral  response  function  for  the  system  was  determined 
using  the  NO  +  0  chemiluminescence  as  an  actinometric 
standard  (Fontijn,  et  al.,  1964?  Keneally,  et  al.,  1967) 
from  4000  A  to  3  g, . 

The  C2Hj  +  0  signal  was  maximized  with  respect 
to  a  set  of  flow  conditions  by  scanning  all  wavelengths 
between  7000  -  9000  A  and  monitoring  the  signal  at  the  contact 
times  determined  by  the  position  of  the  KBr  windows  (contact 
time  in  the  system  at  a  given  pressure  is  from  3  to  15  msec, 
depending  on  the  viewing  port) .  For  this  total  pressure, 
about  1280  g,,  the  first  window  gave  the  largest  signal  at 
all  wavelengths.  The  signals  from  the  other  hydrocarbon  plus 
atomic  oxygen  reactions  were  also  monitored  at  the  first 
window  and  approximately  equal  concentrations  of  hydrocarbon 
(about  20  g,)  and  molecular  oxygen  (about  1200  ^)  were  employed. 
Consequently,  the  intensities  may  not  be  the  maximum  intensitie 
obtainable  although  the  intensities  of  the  signals  from  the 
various  hydrocarbons  correspond  to  the  same  experimental 
conditions. 
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Figure  2.  Intensity  of  chemiluminescence  in  the  700-900  nm  region 
from  0  atom  attack  on  saturated  hydrocarbons. 
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Figure  3. 


Intensity  of  chemiluminescence  in  the  700-900  nm  region 
from  0  atom  attack  on  unsaturated  hydrocarbons. 
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Potto  of  1 


Figure  4.  Intensity  of  chemiluminescence  from  simultaneous  0  atom 
attack  on  and  NO  in  equal  concentrations  of  20y. 
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Results  and  Discussion 


The  chemiluminescent  spectra  generated  by  the  reaction  of 
oxygen  atoms  with  a  series  of  saturated  and  unsaturated 
hydrocarbons  presented  in  Figures  2  and  3  along  with  the 

relative  spectral  response  of  the  photomultiplier  - 
spectrometer:  system,  as  determined  through  calibration  with 
NO  +  0.  At  the  first  window,  the  signal  was  scanned  over  the 
sensitive  wavelength  range  of  the  photomultiplier.  The 
spectra  of  olefin  plus  oxygen  atom  reactions  exhibited  the 
0H(A2£+  -  X2tt)  band  at  3064  A  and  the  emission  in  the  7000- 
9000  A  region  was  characteristic  of  the  OH  Meinel  bands  plus 
a  strong  band  near  7620  A  which  could  be  02(1^)  emission. 

The  intensity  of  the  Meinel  bands  is  a  strong  function  of  the 
contact  time  whereas  the  02(1]£)  emission  is  much  less 
sensitive  over  the  contact  time  range  of  this  experiment. 

The  O  plus  reactive  hydrocarbon  emission  is  structured, 
and  nearly  as  strong  as  the  NO  +  O  emission  in  the  wavelength 
region  from  7000  -  9000  A.  Figure q4  presents  the  spectrum 
of  O  +  C2H4  +  NO  from  7000  -  8000  A.  The  C2H4  +  O  emission 
peaks  are  clearly  distinguishable  superimposed  on  the  NO  +  O 
continuum  emission. 

While  acetylene  exhibits  the  same  structure  in  the 
7000  to  9000  A  region  as  do  olefins,  the  spectrum  from  2500 
7000  A  is  quite  different  with  strong  emission  near  6000  A, 
which  should  be  readily  observable  in  the  presence  of  NO  as 
well,  provided  that  appropriate  spectral  resolution  is 
provided . 

Currently,  we  are  examining  the  reaction  mechanism  of 
jjiki  for  consistency  with  u.v.,  visible  and  infrared  observations 

for  reaction  of  0  +  C2H4: 
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0  +  c2H4 

- 

CH3  +  CHO 

ch3  +  0 

-4 

HCHO  +  H 

HCHO  +  0 

-4 

CHO  +  OH* 

CHO  +  0 

- 

CO  +  OH* 

CHO  +  0 

- 

C02  +  H 

OH  +  0 

- 

0*  +  H 

CHO  +  H 

-4 

H2  +  CO 

The  OH  emission  steps  can  be  initiated  either  using 

directly  or  by  using  HCHO.  Both  generate  strong  2.7ii 

emission.  HCHO,  however,  does  not  seem  to  lead  to  strong 
o 

3064  A  or  Meinel  radiation.  Installation  of  our  Extranuclear 
Quadrupole  Mass  Spectrometer  in  the  next  funding  period  will  greatly 
enhance  our  capability  for  interpreting  these  mechanisms. 

In  order  to  generate  OH  concentrations  suitable  for 

reacting  and  subsequent  observation  of  chemiluminescence, 

H2  is  discharged  using  a  2450  MHz  microwave  system  as  outlined 

in  the  previous  section.  The  concentration  of  TH]  produced 

* 

in  the  discharge  can  be  determined  from  the  H  +  NO  -»  HNO  +  hv 
o 

(7460  A)  titration  as  shown  in  Figure  5.  Normal  [Hi  concentrations 
in  a  total  pressure  of  about  1  torr  are  3  to  8%  (Figure  5  shows 
a  6%  concentration) .  At  the  second  inlet  port  N02  is  added  to 
produce  OH  via  H  +  N02  -  NO  +  OH.  Without  further  addition 
of  reactants,  a  substantial  2.7  and  6.3^  signal  is  observed; 
presumably  this  is  due  to 

OH  +  OH  -  HjO*  +  0 
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FIGURE  9 
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but  some  contribution  may  come  from  OH  +  H2  -*  ^O  +  H. 

Once  the  OH  is  generated,  various  simple  hydrocarbons  can  be 
added  and  the  chomiluminescense  in  the  IR  determined.  To 
date  we  have  done  preliminary  measurements  on  the  added 
reactants  CO,  C2H4  and  C^O  using  IR  filters  and  a  range 
of  solid  state  detectors  (described  earlier) .  Our  results 
show  that  the  OH  +  CH20  reaction  leads  to  a  large  2.7g,  signal 
and  a  somewhat  smaller  signal  in  the  5  -  lOg,  region  (presumably 
HO  bending  mode).  The  OH  +  C2H4  gives  a  rather  small  2.7g, 
signal  but  a  larger  signe  d  in  the  5  -  10  g,  region  than  from  CH20. 
With  CO  +  OH  we  have  not  been  able  to  identify  any  significant 
IR  radiation  although  there  appeared  a  hint  of  a  signal  in 
the  12  -  20 g,  region  (C02  bending  mode)  .  Using  the  available 
actinometric  standard  0  +  NO  at  2.7 g,,  we  have  deduced  a 
preliminary  rate  coefficient  for  the  reaction 

OH  +  CH20  -  2.7gt  photons. 

Assuming  that  the  2.7^  radiation  is  from  ^0*,  we  note  that 
the  system  is  radiaticn  limited  due  to  the  relatively  high 
concentration  of  H2  (present  from  the  discharge  flow)  and  the 
vary  high  vibrational  relaxation  times  associated  with  H20 
(roughly  1  g,sec  at  1  torr  H2)  .  Using  these  numbers,  we 
estimate  that  the  2.7g,  production  from  OH  +  C^O  has  a  rate 
constant  between  2~12  and  5‘14  depending  on  the  [OH]  concentration 
in  the  flow  tube.  We  are  currently  working  on  a  numerical 
model  for  the  flow  tube  chemistry  which  will  enable  a  considerably 
more  reliable  estimate  of  the  [OH]  and  other  intermediate  species 
for  use  in  deducing  rate  constants  and  photon  yields  from 
measured  flow  tube  data. 

These  studies  will  continue  jointly  with  the  development 
of  the  flow  tube  chemistry  code  for  the  interpretation  of  the 
chemiluminescense. 
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Theoretical  Cross  Sections 


In  a  series  of  papers  beginning  in  1964  and  earlier 
Convair  Reports,  R.  Marriott  has  generated  theoretical  cross 
sections  for  many  vibrational  relaxation  processes  involving 
HjO,  C(>2 ,  CO,  Nj,  H2  and  0  atoms.  To  date,  many  of  these 
theoretical  estimates  are  the  only  available  results  to  use 
in  model  non- equilibrium  calculations  involving  these 
molecules.  Yet  these  cross  sections  have  not  in  general 
been  available  in  a  convenient  form  for  use  in  model  cal¬ 
culations. 


To  eliminate  this  problem,  we  have  written  a  simple 
numerical  program  to  reduce  the  cross  sections  calculated 
for  all  vibrational  processes  analyzed  by  the  close  coupling 
method,  including  excited  vibrational  state  processes,  to 
temperature  dependent  rate  coefficients.  The  results  of  this 
analysis  will  be  described  in  a  forthcoming  report  where 
available  experimental  data  will  be  included. 


The  results  of  this  analysis,  particularly  the  excited 
level  rate  coefficients,  shouldbe  important  in  providing 
scaling  information  of  the  relaxation  of  excited  vibrational 
levels  and  also  provide  currently  unknown  relaxation  rates  for 
HjO  (stretching  modes  and  bending  modes)  in  collision  with  H2, 
C02  and  H20  as  collision  partners.  The  relaxation  processes 
being  calculated  are  tabulated  in  the  following  list. 


1?  C02(v2)  j 

C02  ( Vl}  ( 

C02 (2v2) J  + 

C02(v3)  ) 


H2 

H2 

CO 


0 

2 


2?  H20(v2)  j  H2 

H20(2v2M  +  H^O 

H2°<V 

H20(V3)  ) 


3.  CO  +  H20,  0(1D),  He 


4.  N2  +  AT 


♦plus  all  interstate  combination  processes 
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SECTION  III 


Molecular  Beam  Studies 

An  experimental  effort  is  being  initiated  to  provide  input 
data  needed  for  calculating  excitation  cross  sections.  Of  partic 
ular  interest  are  those  with  atomic  oxygen  at  energies  between 
a  few  tenths  of  an  electron  volt  to  several  electron  volts. 

An  apparatus  designed  for  this  type  of  experiment,  and 
constructed  by  High  Voltage  Engineering  Corporated  for  NASA 
Lewis  Research  Center  was  donated  to  Wayne  State  University 
when  the  NASA  program  was  discontinued.  It  is  proposed  to 
use  this  equipment  to  measure  the  needed  potentials. 

The  equipment  has  been  transferred  to  our  laboratory, 
and  during  the  present  reporting  period  work  has  been  carried 
out  to  bring  the  machine  back  to  a  functioning  level,  and  to 
install  the  necessary  modifications  for  these  measurements. 

The  machine  is  now  almost  ready  for  preliminary  experiments. 
A  post-doctoral  fellow  will  join  our  group  in  the  middle  of 
January  to  assist  in  carrying  out  the  research  program. 
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